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MESOSCOPIC MODELING OF STOCHASTIC REACTION-DIFFUSION KINETICS IN 

THE SUBDIFFUSIVE REGIME 


EMILIE BLANC, STEFAN ENGBLOM, ANDREAS HELLANDER, PER LOTSTEDT* 

Abstract. Subdiffusion has been proposed as an explanation of various kinetic phenomena inside living cells. In order to 
fascilitate large-scale computational studies of subdiffusive chemical processes, we extend a recently suggested mesoscopic model 
of subdiffusion into an accurate and consistent reaction-subdiffusion computational framework. Two different possible models 
of chemical reaction are revealed and some basic dynamic properties are derived. In certain cases those mesoscopic models 
have a direct interpretation at the macroscopic level as fractional partial differential equations in a bounded time interval. 
Through analysis and numerical experiments we estimate the macroscopic effects of reactions under subdiffusive mixing. The 
models display properties observed also in experiments: for a short time interval the behavior of the diffusion and the reaction 
is ordinary, in an intermediate interval the behavior is anomalous, and at long times the behavior is ordinary again. 

Key words, continuous-time random walk, subdiffusion, fractional derivative, anomalous kinetics, multistate reaction- 
diffusion system. 
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1. Introduction. Quantitative models of reaction-diffusion systems are important tools to theoretically 
and computationally study the dynamics of intracellular control systems. Macromolecules such as proteins, 
mRNA and DNA interact in regulatory pathways in which each chemical reaction occurs in a certain part of 
the cell. Molecules are transported via diffusion or active transport to arrive at the subcellular location needed 
to fulfill their function. As an example, in many gene regulatory pathways, proteins called transcription 
factors will diffuse and search for the correct binding site on DNA, where they modulate the translation 
of DNA into mRNA, which in turn will be translated into proteins. Apart from the spatial dynamics, 
realistic models need to account for large fluctuations in the copy numbers of the species due to the small 
reaction volume and thus the small number of molecules of key species such as transcription factors. To 
that end, spatial stochastic models based on a Markov process formalism are popular due to their high level 
of biological realism compared to macroscopic partial differential equations (PDE), with only a moderate 
increase in computational complexity and cost. Such mesoscopic models, based on the so called Reaction- 
Diffusion Master Equation (RDME), have recently been used to address different biological phenomena such 
as regulation of cell division in E. coli [11], yeast polarization [31] and genetic oscillators [50]. 

A fundamental assumption in the RDME models is that molecules are point particles. Thus, the model 
does not account for molecular volume exclusion, something that can lead to anomalous diffusion in the 
crowded compartments of living cells. Especially on 2D membranes, crowding and anomalous diffusion 
behavior can be expected to have an impact on the reaction kinetics. For diffusion limited reactions, sub¬ 
diffusion often results in a slower decay [54]. A reversible ligand binding bimolecular reaction on a 2D 
membrane is studied in [46] where anomalous diffusion is simulated with a continuous-time random walk 
(CTRW). The steady state distribution of the bound complex depends critically on the anomalous diffusion 
and its parameters. An irreversible bimolecular reaction is simulated in [1] with ordinary diffusion and a 
modihed non-Fickian diffusion which is anomalous for intermediate time but tends to ordinary diffusion in 
the long term. The distribution of the reactant is sensitive to the type of diffusion. The conclusion in [40] 
is also that subdiffusion changes the reaction rate in a bimolecular reaction. 

To capture effects of subdiffusion due to molecular crowding in stochastic reaction diffusion simulations, 
detailed Brownian Dynamics (BD) models based on a hard-sphere assumption can be used instead of the 
mesoscopic RDME model, since they account for the volume exclusion of molecules. However, this comes 
at the price of a very large increase in computational cost, making simulations on the timescales of interest 
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in systems biology (minutes to hours) challenging. Thus there is a need to investigate less computationally 
demanding approximations of mesoscopic subdiffusion. A promising approach for mesoscopic anomalous 
diffusion simulation was suggested recently by Mommer and Lebiedz [38], in which a stochastic CTRW 
model [39] is approximated by an internal states model where all transitions between states have exponentially 
distributed waiting times. Hence, the subdiffusive process is simulated as a coupled reaction-diffusion system 
with ordinary diffusion, making it easily implemented in software frameworks based on the RDME [8, 18, 22]. 

Starting from the viewpoint that the mesoscopic model arises as an approximation of molecular crowding 
on the microscopic level, the internal states model can be viewed as a mathematical means to arrive at an 
approximate, faster simulation method. On a more general level, such internal states models arise also from 
the need to model for example conformational changes of macromolecules. Many proteins exist in a number 
of different states due to e.g. ligand binding, methylation, and conformational changes [4, 9]. In practical 
modeling, if each state is to be represented as a separate chemical species, an exponential growth in model 
complexity due to a combinatorial explosion in the number of states may be the consequence. To handle 
this from a practical point of view, special rule-based modeling languages have been developed, such as 
BioNetGen Language [3] and PySB [33]. This combinatorial explosion in states motivates the development 
of computational methods to simulate models with these types of molecules more efficiently [49]. In this case, 
we start with an internal states model naturally arising from the biological model, and are then interested in 
how this model can be approximated on a phenomenological, macroscopic level as a fractional PDE (FPDE). 

Molecules transported by subdiffusion move slower than with ordinary diffusion [36] which could be an 
effect of crowding in a biological cell. The governing equation for subdiffusion is a FPDE with a fractional 
time derivative and the Laplacian as the diffusive space operator. Chemical reactions can also be included 
in this framework. Three different models for monomolecular reactions with subdiffusion are derived and 
compared with analytical solutions in [21]. In one of the models where the fractional derivative acts only on 
the diffusion, the concentration of the species can become negative, making it less suitable. In the preferred 
model, where a special fractional derivative is derived for the diffusion, the equation for the homogeneous 
solution is recovered without a fractional derivative. The model where the fractional derivative is applied to 
both the diffusion and the reaction is used in [53] for a bimolecular reaction. A model with internal states is 
derived in [43] for the propagation of a reaction front in an inhomogeneous medium. The particles move by 
diffusion and react with each other at the front. Macroscopic equations are obtained after summation over 
the internal states. Analytical results show that the speed of the front depends on the type of reaction and 
the particle distribution between the internal states. That same model is suggested in [44] to explain the 
results in [46]. 

In this paper, we develop theory to extend the mesoscopic anomalous diffusion internal states model 
[38] to account also for chemical reactions, making it possible to model general reaction-diffusion processes. 
In particular, we develop theory for the connection between the model framework that thus emerges and 
existing models of anomalous reaction-diffusion processes on the macroscopic, FPDE level. As we will see, 
the proposed reactive internal states model is quite general and can, depending on the involved parameters, 
result in different mean-field equations. Some of these mesoscopic equations have a simple interpretation 
at the macroscopic level as a FPDE but in general this is not possible. In numerical simulations, the 
reaction-diffusion systems have the same behavior as observed in many experiments: in a short time interval 
after start we see ordinary behavior, then there is an anomalous phase, and finally the behavior is ordinary 
again but with different diffusion and reaction coefficients. These conclusions are supported theoretically 
in [23, 25, 30, 41]. The reaction-diffusion equations encompass many of the known models of subdiffusive 
reaction systems, making it promising as a general modeling framework. 

The rest of the paper is organized as follows. In Section 2, we review some basic properties of stochastic 
models of diffusion and how they give rise to a phenomenological, macroscopic FPDE in the thermodynamic 
limit. Here, we will also see how the addition of chemical reactions on the macroscopic level can lead to two 
different, possible FPDE models. We introduce an internal states model of reaction-subdiffusion systems in 
Section 3 and analyze its mean field properties. In Section 4, we relate the models in Sections 2 and 3 to each 
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other for linear and non-linear chemical reactions, and in Setion 5, numerical experiments are presented. 
Finally, the paper is concluded in Section 6. 

2. Fractional partial diflferential equations as limits of continuous-time random walks. In 

this section, we recall some well-known facts about Brownian and, in particular, subdiffusive random move¬ 
ments. When approached in the proper macroscopic limit, FPDEs emerge as a convenient mathematical 
model for subdiffusion. Terms are added to the FPDEs to model chemical reactions. 

The CTRW model was introduced by Montroll and Weiss for hopping transport on a disordered lattice 
[39]. In this model in continuous space, the particle is assumed to traverse the space by a series of jumps. 
The displacement and the waiting time to perform the next jump are drawn from a given probability 
density function (PDF) We assume that the jump length PDF X{x) and waiting time PDF are 

independent random variables. Consequently, is written 

'^{x,t) =tp(t)X{x). (2.1) 


Different diffusion processes can be categorized by the expected waiting time 


and the jump length variance 


dt, 





\\x\\lX{x)dx, 


( 2 . 2 ) 


(2.3) 


where d is the dimension of the embedding space. If both and r* are finite, the long-time limit corresponds 
to Brownian motion. A diverging r* with hnite gives rise to subdiffusion. On the contrary, a diverging 
E^ with a finite t* induces superdiffusion [23, 36], which is beyond the scope of this article. 

2.1. Brownian motion. We consider a Gaussian jump length PDF 


A(a;) 


(4 7ra2)‘^/2 


(2.4) 


and a Poissonian waiting time PDF 

^(t) = i 

r 


(2.5) 


Equations (2.2)-(2.3) lead to E^ = 2 cr^ < oo and r* = r < oo. Since E^ and r* are finite, the long-time limit 
thus corresponds to Brownian motion. At the macroscopic scale, we recover the classical diffusion equation 
[23, 36] for the concentration U of the chemical species A 


with 


f—■ 



The linear time dependence of the mean squared displacement 


{\\x\\l{t)) = 2dDt 


( 2 . 6 ) 


(2.7) 


( 2 . 8 ) 


is characteristic of Brownian motion. 
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2.2. SubdifFusion. We consider now instead a Gaussian jump length PDF (2.4) and a waiting time 
PDF [28, 35] 


where 


V'(t) = 




■E^ 


0 < a < 1, 


= E 


i-jt/rrr 

F(a k + a) 


(2.9) 


( 2 . 10 ) 


is the generalized Mittag-Leffler function. The Mittag-Leffler waiting time PDF has been observed ex¬ 
perimentally, from polymer rheology [16], over ligand rebinding to proteins [15] and protein conformation 
dynamics [52], to financial market time series [34]. Equations (2.2)-(2.3) lead to = 2 < oo and r* = oo, 

which give rise to subdiffusion. At the macroscopic scale, we obtain the FPDE [23, 36] 


d°‘U 

dt°‘ 


= AC/, 


with 


if - E - 

1\ rv - - 


2r° 


( 2 . 11 ) 


( 2 . 12 ) 


At the boundary dfl of the domain D, the molecules are reflected back implying homogeneous Neumann 
boundary conditions at dfl. The operator involved in (2.11) is a Caputo fractional derivative in time of 
order a, generalizing the usual derivative. It is defined as [5, 37] 


5“t/ 

dt°‘ 


{t - t)-“ dU 

r(l — a) dt 


(r) dr. 


The mean squared displacement is given by the power law 


(11=1=112(0) = 


2dKa, 


■ r. 


(2.13) 


(2.14) 


r(l -f a) 

In the one dimensional (ID) case (d = I) in free space, the Green’s function of the FPDE (2.11) is [36] 


U{x, t) = 


•^4 TT Ka 


rr2.0 

^1,2 



(’-f-) 


4a:„*“ 

(0,1) 

-1 

1—1 


(2.15) 


where the Fox function is defined in Appendix A. In the particular case a = 1/2, the Green’s function 
(2.15) can be rewritten 


U{x,t) = 




^3.0 / X 

I 16 7 ^ 1 / 2 * 1/2 


1 1 


(2.16) 


where the Meijer-G function Gq’j is also defined in Appendix A. The Green’s function (2.16) will be used in 
numerical experiments in Section 5 as a reference solution. 

2.3. Adding chemical reactions. In the following we discuss how chemical reactions can be added 
to the macroscopic subdiffusion model (2.II) for the cases of annihilation, reversible isomerization, and 
reversible bimolecular reactions. As we will see, in the first two cases two different FPDE models are 
possible, while only one of them turns out to be well defined for bimolecular association. 
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2.3.1. Annihilation. We first consider one species A and the annihilation process A — A0. Two dif¬ 
ferent FPDEs can be used to model an annihilation process. If a constant proportion of walkers are removed 
instantaneously at the start of each step then the long-time asymptotic limit yields a fractional reaction- 
diffusion equation with a fractional order temporal derivative operating both on the standard diffusion term 
and on the linear reaction kinetics term [21, 24, 42]: 


dU 9^““ 


(2.17) 


In what follows, (2.17) is referred to as model I. The total amount [/ of A in a bounded domain 17 with 
boundary 317 is defined as 


U = / U{x,t) dft. 


(2.18) 


The equation (2.17) is integrated over 17 

du 31 -“ 


dt 371 -“ 


Kr. 


MJ dVt-KU] . 


Since 


(2.19) 


Af7dl7 = 


n-VUdS = 0 


I on 


( 2 . 20 ) 


for Neumann conditions at 317 with normal n, U satisfies the fractional ordinary differential equation (ODE) 

d°‘U 


dt° 


= -K U. 


( 2 . 21 ) 


Consequently, the total amount of A is also affected by the subdiffusion. This behavior is called anomalous 
kinetics. In the ID case in free space, the Green’s function of (2.17) is [21] 


1 {—k^t°‘y 2,0 


(l-|+aj,a) 


V47rA„7« j! 

1 

(0,1) 

(^+Al) . 


U{x,t) = 

In the particular case of a = 1/2, (2.22) can be rewritten as 

1 


( 2 . 22 ) 


U{x,t) = 


j=o 


(—2 k^, t^l'^y ^3 0 


j! 


0,3 


16Ai/2 71/2 


1 I 1 

4 2’ 2 


(2.23) 


If instead the walkers are removed at a constant per capita rate during the waiting time between steps 
then the long time asymptotic limit has a standard linear reaction kinetics term but a fractional order 
temporal derivative operating on a nonstandard diffusion term [21, 32, 45, 51]: 


— (e'^* * AC/) - fc, U. 


dt “ 3/1 

This case is referred to as model II in what follows. Integrating (2.24) and using (2.20) leads to 


^ -K e 
dt " 


— fc* t 


31 


371 - 


t 


AUdn] - k^U = -k^U. 


(2.24) 


(2.25) 
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Ordinary kinetics is recovered: the total amount of A is not affected by the subdiffusion. The change of 
variables C/ = in (2.24) leads to the FPDE (2.11). In the one dimensional case, the Green’s solution 

of (2.24) is then [21, 36] 


l/(x, t) 


-fc* t 


H 


2,0 


•\/4 TT Ka 



(‘-I’”) 


-1 

1—1 

1—1 

o' 


(2.26) 


2.3.2. Monomolecular reactions. Let two species A and B undergo the monomolecular reversible 
isomerization reaction A^B. We assume that the jump length variance and the time scale involved in 

(2.12) are the same for each species: and ta = tb = t. Consequently, the diffusion coefficients 

are also identical Kau = Ksa = Ka- The concentrations of A and B are U and V, respectively. As in the 
case of annihilation, two different FPDEs can be used to model this reaction. Model 1 is here written 


dU 3 ^““ 

= K^ —— (At/ -RU) 


(2.27) 


with U = {U Vy and R = 


p = 


1 I 

— 1 k 


. The reaction matrix R is diagonalizable, R ^ P Ar P ^, with 

(2.28) 


Ar = 


k^ 0 
0 0 


By setting tj = P ^ U, the governing system of evolution equations decouples to 

dil d^~°‘ t ~ ~\ 

In the ID case without boundaries, the Green’s function of the FPDE (2.29) can be computed using (2.15) 
and (2.22). 

Model 11 is written 

^ AC/) - iCC/, (2.30) 

The change of variables tj = in (2.30) leads to the FPDE (2.11). The Green’s function of the FPDE 

(2.30) in ID can be computed using (2.26). 

2.3.3. Bimolecular reactions. Finally, we consider three species A, B and C and the reversible 

fe. . 

bimolecular reaction A -\- B . As previously, we assume that the jump length variance and the time 

scale involved in (2.12) are the same for all species, which implies that the diffusion coefficients are also 
identical. The concentrations of A, B and C are C/, V and W, respectively. In the bimolecular case, due to 
the nonlinearity, only model 1 is defined. It is written 

dU 9^““ 

= Rc.g^(AU-k^UV + £.W), 

< — = K^g^(AV-k,UV + £,W), (2.31) 

dW 9^““ 

^ = K^Q^{AW + k.UV-l.W). 

No analytical solution is available in this case. Nevertheless, the quantity U — V satisfies the FPDE (2.11). 
Gonsequently, the analytical solution of C/ — P is given by (2.15). 
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3. The internal states model. In this section, we introduce an internal states model of subdiffusion 
on the mesoscale. The molecules jump on a lattice and can react with each other when they are in the same 
lattice cell. Diffusion is here modeled by discrete jumps with a waiting time PDF tjj(t) and a jump length 
depending on the nearest neighbors in the lattice. 

3.1. Mesoscopic diffusion and reactions. In the stochastic mesoscopic model, the state of the 
system at time t is given by the number of molecules of each one of the N species in the state vector 
y{t) G 1^. The domain is partitioned into non-overlapping voxels or compartments Vi, i = 1,... ,M, with 
a mesh or a lattice. Each voxel has a vertex at Xi in the center. The component yjk of is the number 
of molecules of species j in Vk- The states are changed randomly at discrete time points depending on 
chemical reactions or a jump to a neighboring voxel due to diffusion. Between ti and t^+i, y is constant. The 
PDF is p{y,t) for the system to be in state y at time t, and satisfies a master equation [10, 26], commonly 
referred to as the Reaction-Diffusion Master Equation (RDME) in the context of chemical kinetics. 

The jump coefficient Xji for a diffusive jump of a molecule A from Vj to Vi is determined from the 
coefficients of a discretized diffusion equation in [10]. The rate for a jump is \jiyj where yj is the copy 
number of A in Vj. After the diffusive jump, yj and yi are updated: yj := yj — 1, yi := yi + 1. 

A bimolecular reaction between the molecular species A and B produces a C molecule in a voxel in 

A + B^C. (3.1) 

The propensity of the reaction is k times the copy numbers y and z oi A and B. If the reaction takes place at 
ti, then the state before the reaction y(t7) in fh® voxel is changed to y(t^) = y{t~^)+i'r immediately after the 
reaction r, where Uj. is the stoichiometric vector associated with the reaction. In (3.1), y := y — 1, z := z—1, 
and the copy number of C increases by one. 

A realization of the chemical system is generated by the stochastic simulation algorithm (SSA) [13] in 
a Monte Carlo method. In the limit of large copy numbers, the mean values of the concentrations of the 
species converge to the solution of the reaction rate equations [29]. These equations are stated in Section 3.3. 

In the next section, a master equation is derived for the PDF of a system with internal states. Internal 
states of the species are introduced in [38] to model subdiffusion. We will extend their model to include 
reactions between the species. 

3.2. A generalized master equation. As a generalization of the classical chemical master equation 
(CME), we let a chemical species A have a number N of internal states Ai, I = 1,... ,N. The state Ai 
can model a conformational state that can be difficult to observe and that is more or less hidden from a 
practical point of view as in a hidden Markov model. On the molecular level, the hidden states can model 
different things, for example, they can represent different geometrical configurations of the molecule or the 
macromolecule can have different small molecules attached to it at different locations [4, 9]. 

We now consider a general event q which could be a diffusive jump, a change of internal states or a 
reaction. Assume that the molecule in the state j has a PDF for the waiting time with density ipqjit) for 
the event q and that the molecule, due to the event, then changes state to i with probability iTqij such that 
"Y^iT^qij = 1. Although formally ijjqjit) here is slightly different from the one in Section 2, the notation is 
retained since the observed effects are the same. Usually, the waiting time is assumed to be exponentially 
distributed and dependent on reaction rates and diffusion propensities. 

The Laplace transform of a function f(t) is denoted by /(s). An auxiliary function 4>qj{t) is defined by 
its Laplace transform involving the Laplace transform of the waiting time PDF '4>qj{s) 

^qj[s) = s^qj{s)/{I - ^qj{s)) (3.2) 

Gillespie shows in [14] (see also [19, 27]) that there is a generalized master equation for the PDF pi of state 
i 

- T^qji (3.3) 
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If ipqiit) = a^i exp(—agif), then 4>qi{t) = aqiS(t) where S{t) is the Dirac delta and Gqi is the reaction 
propensity. Then we obtain the usual master equation 

dp'(t') 

= 'y ] y ^qj Pji^) ~ '^qji ^qiPiit)) = y ] 'y^i'^qij 0>qjPj{ty) — aqiPi{t). (3.4) 

q j q 3 

Let state * be y and let state j be y — Uq in (3.4) and identify ag(y — Uq) with iTqij Uqj and Uqi with 
o-qiy) for event q. Let Q be the total number of events. If the waiting times are exponentially distributed 
then we recover the usual CME [26] 

= '^iaq{y - iyq)p{y - i^q,t) - aq{y)p{y,t)). (3.5) 

9=1 



event 

probability 

propensity 

Change of internal state 

Diffusion from voxel k to voxel 1 
Monomolecular reaction 
Bimolecular reaction 

Production of one molecule 
Annihilation 

Aj — ^ Ai 

Ajk Ajj 

Aj — >■ Bi 

Ai + Bj — >■ Ck 

0 ^ Ai 

Ai ^ 0 

T^ik = ^ik/ Afc 

TTij = Kji/ Kj 
Itij — Kijk / Kij 
'^iO — 

T^Oi = 1 

oj = %7d 

Afc Vjk/Tj 

Vf 

Qij — Kij Pi Zj jTij 

ai = kjn 
ai = kpj/n 


Table 3.1 

Jump probabilities and propensities in (3.4) for different events. 


We consider six different elementary events listed in Table 3.1. The coefficients in (3.4) and waiting 
times in (3.2) defining the events are also found in Table 3.1. The copy numbers of the species Aj and Bi 
are pj and Zj, respectively. The total rates away from the state in the table to all other states that can be 
reached are 

y ( pi — 1? — y ( ^tk: tvj — y ( ktjij kiij — y ( tvijk 

i I i k 

ensuring that = 1- Th® waiting time distribution in the state is ip jit) = aj exp(—Oj t) when one 

molecule is transformed in Table 3.1. The time scale Tj of the transformations depends on the internal state 
j of A. We will find that if ^ Tj when i ^ j then the annihilation, diffusion, reaction, or production is 
anomalous at the macroscopic level. If = t then the transformation between the states is the ordinary 
one with the same waiting time for all states. 

The scale of the bimolecular reaction depends on the states of A and B. By splitting the reaction 
rate into OKijk/Ti + (1 — d)K,ijk/Tj with 0 ^ 0 ^ 1, the time constant is = l/(0r“^ + (1 — 9)Tp). This 
Tij is in agreement with the time scale obtained from [6]. Let DAi and Dsi be the diffusion coefficients of 
species A and B in state i, kb the microscopic reaction rate, and pr the reaction radius. Then the reaction 
rate Kijk in [6] follows from Smoluchowski’s rate law and is a function of Dai and Dsj 

_ kb4:TriDAi + Osj) Pr 

kb+AT{DM+DB,)pr' ^ ’ 

Let the diffusion coefficients of species A and B depend on the internal state such that Dai = cr\lTi and 
Dsj = (^%lTj (2.7). Then with 6 = <j\/ia \-\-and for diffusion limited systems with a large kb compared 
to DAi + Dsj, Kijk is approximated by 


Kijk ~ 4 TT (ct^ + a'g) 


e 1-6 

-1- 

n Tj 


(3.7) 











When kh is small then the influence of the diffusion disappears in (3.6) and Kijk ~ fcfe without dependence 
on the states i and j. 

With the master equation (3.4) or (3.5), we can derive reaction rate equations approximately satisfied 
by the mean values as in [26, 38]. 

3.3. Mean-field properties of the internal states model. The mean values of the concentrations 
of the species approximately satisfy a system of ODEs often denoted the reaction rate equations. These 
equations are derived from the PDF in the master equation (3.4) or (3.5), see [26]. When the reactions are 
such that all propensities are linear in the chemical system, see Table 3.1, then the solutions to the equations 
are the exact mean values. They are approximations if there is a bimolecular reaction in the chemical system. 

3.3.1. Diffusion. Let us first examine a system with one molecular species and N internal states. This 
is the problem investigated in [38]. We allow changes of internal state and diffusion between two voxels but 
without chemical reactions. The mean values of the concentrations of the N states Ui{t) G in voxel i are 
the solution of 




A\ii. 


(3.8) 


The number of voxels directly connected to Vi is rii. Consequently, a diffusive jump between Vj and Vi is 
possible. The elements of A follow from Table 3.1 and are Aij = fJ-i/rj, i ^ j, and An = (/x^ — l)/ri. Let T 
be a diagonal matrix with Tn = I/tj, p a vector with non-negative components pi such that e^p = 1 and 
e a vector with Cj = 1 for all i. Then A in (3.8) can be written 

A = - I)T. (3.9) 

The nullspace consists of one vector Uqo where 

Uoo=T~^H. (3.10) 

The left eigenvector of A corresponding to eigenvalue Ai = 0 is e such that 

e^A = 0. (3.11) 

The diffusion jump coefficients in (3.8) are derived such that the Laplacian is approximated in voxel i 


Au(a;, f) ~ ^ Ay - A^ u^. 
1=1 


(3.12) 


On a Cartesian mesh with equal mesh spacing h, Xji = 1/h^ and Xi = 2d/h'^ where d is the dimension. 
On an unstructured mesh, the coefficients can be derived by a finite element method as in [10]. With a 
continuous u(a;,t) in space, the equation approximated by (3.8) is 


/^ll 

— =DAu + Au, D = a^T. 

at 


(3.13) 


The boundary conditions are of Neumann type at the boundary to preserve the total concentration of the 
species as in (2.20). The analysis is simplified in this section if we consider the solution u(a;,t) of (3.13) 
instead of the discrete solutions Ui(t) in the voxels in (3.8). 

With positive initial data, it is easy to see that there is a unique positive steady-state solution Uqo as 
t —1 oo. This solution is space independent. We write this in the normalized form 


lim u{x,t) = Uoo = uuoo, 

t—VOO 


(3.14) 


9 


where ||uoo||i = 1 and u = ||u(a;,0)||i by the preservation of mass. Since Auoo = 0 we have Au^o = 0 and 
Uoo is given by (3.10) and (3.14). 

Expand the solution of (3.13) in a cosine series in ID 


OO 

VL{x,t) = ^U„(t) cos(wx), 


with X in [0, 27r] and insert into (3.13). Then each mode satisfies 


The solution to the equation is 


u^(t) = S exp(At) S ^ Uui{t), 


(3.15) 


(3.16) 


(3.17) 


where S{oj) = (si, S 2 ,..., s^r) is the eigenvector matrix of —uj^ a^T + A and the eigenvalues Xj{u}) are on the 
diagonal of A. By Gerschgorin’s theorem for the eigenvalues of a matrix, the eigenvalues all satisfy 5iAj ^ 0 
for w = 0 and 5iAj < 0 for w > 0. Thus, as t increases all modes vanish except for one mode si(0) = Uoo at 
w = 0 with eigenvalue Ai = 0. 

The concentrations of the internal states are summed at the macroscopic level. Then U = e^u satisfies 


dU{x, t) 

m 


= e'^{a^T Au + A u) = T Au = cr^ T (t) cos(wa;) 

= ^{x, t) Au = j{x, t) AU, 


(3.18) 


where 


^ T (t) cos(a;a;) 

j{x, t) = ^ -. (3.19) 

^ u,^(t) cos(a;x) 


The macroscopic U satisfies a diffusion equation with a diffusion coefficient varying in space and time. For 
large t, the dominant mode in the spatially non-constant part of the solution is damped by the eigenvalue 
Ai = maxAj(a;) < 0 at a;i. Then Si(a;i) and the steady-state macroscopic diffusion coefficient in 

(3.19) is 


2 uJi T Si cos{iUix) 22^1 

7 ~ rr — 7^ — 7^ -^^^— = (j —— -= (j 


N 

sii/n 


e'^ Si 


ujf Si cos(wix) 

3.3.2. A reversible reaction. Consider next the simple case of a reversible isomerization, 


N 

E sii 

i=l 


(3.20) 


A,; 


(3.21) 


which is to be understood in the sense that, for example, the rate for the jth internal state of a i?-molecule 
to transform into the ith state of an A-molecule is Xji. These are thus two monomolecular reactions with 
reaction rates in Table 3.1. 
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We assume that the variance of the jump length, the diagonal matrix T and the vector jj, are the same 
for all species with identical matrices D and A. Taking the mean as in (3.8) we readily arrive at the coupled 
set of PDEs 


— = D Au + A u — Ki u + L 2 V, (3.22) 

9v 

— = D Av + A V + K 2 u — V, (3.23) 

for some matrices Ki, K 2 , and Ti, L 2 whose precise form we now determine. Define positive rate matrices 
K and L with Kij = Kij and Ly = Ay. From the prescription (3.21) we find that in (3.22)-(3.23), the ith 
state is affected by the reaction terms 


N N 

( ~t“ L 2 v)2 — ^ ^ y^i T ^ ^ Aji , 

i=i 

N N 

(~t“.K^2 Li v )2 — -j- ^ ( Kji llj ^ ( Ay' Vi- 
3 = 1 3 = 1 

Identifying terms we obtain that 

Ki = diag(li:e), K 2 = K^, 

Ti = diag(Xe), L 2 = , 

where the notation diag(f) denotes a diagonal matrix with fi on the diagonal. 

We are interested in the stable, space independent solutions to (3.22)-(3.23). Introduce 

A-Ki L2 

K2 A-Li 

By construction and using (3.11) we have the crucial properties that 

M 

k = l 
k^i 

Bij ^0, j. 



(3.24) 

(3.25) 


(3.26) 

(3.27) 


(3.28) 


(3.29) 

(3.30) 


with M = 2N. Also, by inspection B is irreducible, i.e. there is no permutation matrix Q such that 

QBQ-^ = (^ Q ^ ^ . (3.31) 

Taken together, B is a W-matrix and this provides us with certain general stability properties. 

Lemma 3.1 {W-matrix lemma). Suppose that a real, irreducible matrix H € ]^MxM gdfigjlgg (3.29)- 
(3.31) (with B replaced by H). Then the system of ODEs 

i'{t) = Hm (3.32) 

has a single stable equilibrium solution ^00 ^5 t —> 00 . Moreover, if initial data ^(0) with positive mass 
e^4(0) > 0 JS given, then e^^{t) = e^^{0) > 0 for all t > 0. 
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This particular formulation is discussed in detail in [26, Chap. V.3] and we note that it can also be shown 
to follow from the Perron-Frobenius theorem. To get some further insight into the stability we present some 
alternative arguments as outlined in [26, Chap. V.9, p. 129]. 

Proof. Consider the adjoint problem 

= (3.33) 


such that 




(3.34) 


Using (3.29) we see that 


M 

fc=i 

kjii 

By the irreducibility of H, the largest element in rj decreases and the smallest element increases with time 
such that an all-constant vector rjoo emerges in the limit. With rjifS) = e^, the ith unit vector, we recover 
from (3.34) the unique equilibrium solution ^oo- Also, with r7(0) = e we have r](t) = e and mass is a 
preserved quantity. □ 

Apply the lemma to B in (3.28) with = (u^,v^). Then it follows that there is a steady state 
(uoo , Voo) and the initial mass e^u(O) -f e’^v(O) is conserved for all t. 

We now turn our attention to the equivalent reaction rates as induced by the subdiffusive reactions 
(3.21). We write the unique equilibrium solution UooWoo in the normalized form Uqo = uUqo, Voo = Woo, 
where ||uoo||i = ||voo||i = 1. Using (3.26)-(3.27) we readily find in good agreement that 

fee, :=e'^/CiUoo =e'^K2Uoo, := ii Voo = ^2 Voo- (3.36) 

Also, to mention just one example, 

fceqe^Uoo = e^iTiUoo, (3.37) 

which can be understood as a kind of consistency result; at steady state the equivalent reaction rate applied 
to the sum of the internal states gives the same result as the sum of the individual subdiffusive rates. 

To conclude, using the subdiffusion steady state solutions, which remain valid also when the coupling 
reactions (3.21) are ‘turned on’, we can read oft the equivalent reaction rates as a function of the corresponding 
subdiffusive rates. Insert u(a;, t) = u{x, t) Uoo and v(a;, t) = v{x, t) Voo into (3.22) and (3.23). Then we have 

£) Uoo Am-I- u Auoo — Uoo + ^2 Voo, (3.38) 

D Voo Am- b M A Voo + MiT2 Uoo - uTi Voo. (3.39) 

Multiply by and in the long-time limit we thus arrive at the familiarly looking macroscopic reaction- 
diffusion PDF 

dvL 

— = -fuAu - keqU + leqV, (3.40) 

dv 

— =JyAv + keqU- leq M, (3-41) 

with 'fu = D Uoo and 7 „ = Voo . 
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du 

'di 

dv 


3.3.3. A bimolecular reaction. We finally consider the case of a reversible dimerization, 


A, + Bj"^"Ck. (3.42) 

^ijk 

As in Section 3.3.2, the jump length variance, the diagonal matrix T and the vector fi are the same for each 
species and D and A are also identical. Due to the nonlinearities it is inconvenient to write this in a matrix 
form as in (3.22)-(3.23). Instead, for each internal state i we have the mean-field equations. 


dui 

~dt 

N 

= Di Am* -f ^ Ay Uj 

N 

^ ^ ^ijk '^i '^j H" 

N 

^ ^ Lijk UJkt 

(3.43) 

i=i 

j,k^l 

j,k^l 


dv^ 

dt 

N 

N 

N 


= Di Avi -f ^ Aij Vj - 

- ^ Kjik Uj v^ + 

^ ^ LjikUJki 

(3.44) 


j,k^l 



dwi 

N 

N 

N 


= Di Awi + ^ Aij Wj 

“h ^ ^ Kjki Uj Vk 

^ ^ LjkiUJi^ 

(3.45) 

i=i 





where for readability we write Kijk = Kijk and Lijk = Xijk- The difference Ui — Vi satisfies 


d{ui - Vi) 


N N 

DiA{ui - Ui) + ^ Ajj (ui - Vi) - ^ KijkUiVj - Kj^kUjVi. 


(3.46) 


If Kijk = Kjik and Lijk = Ljik in (3.46) and Ui{x, 0) = Vi{x, 0) for all i, then it follows that Ui{x, t) = Vi{x, t) 
for f > 0. 

Assume that there is a positive steady state solution (uqo ,Voo ,Woo) to (3.43)-(3.45) when t ^ oo. Let 
(Uoo, Voo, Woo) = (uUoo ,uVoo , wwoo) where |luoo||i = ||voo||i = ||woo||i = 1. Small perturbations around 
the equilibrium solution are evolved by the Jacobian of the system. For (3.43)-(3.45), it can be written in 
the form 


A - ATn -Ki 2 Li 
B = I -K 21 A - K 22 L 2 
K 31 K 32 A — L 3 


We identify after some tedious work 

N 


N 


N 


k=l 


k=l 


(3.47) 


Kii = diag ^ Kijk 

Uooj-i K .22 — diag ^ ^ LjikUooj^ 

L 3 = diag ^ Ljki, 

(3.48) 

j,k^l 

3 ,k^l 

j,k=l 


N 

N 



(-^ 1 ) 2 ^ = Ljkj , {L 2 

)u “ Lkij , 


(3.49) 

k^l 

k^l 



N 

N 



(-^12)zj - ^ ^ Kijk Uoci-; 

{K 21 ^ij — ^ ^ Kjik Uoci^ 


(3.50) 

k^l 

k^l 



N 

N 



{K^ 2 ')ij — ^ ^ Lkji Uocki 

(-^ 31 ) 2 ^' — ^ ^ Kjki Uook- 


(3.51) 


The Jacobian B in (3.47) is not a W-matrix. To understand why, note in (3.42) that the total sum of 
molecules is not a preserved quantity. However, if each C-molecule is counted twice, this new weighted sum 
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is in fact preserved. Indeed, by inspection and after some work we find that multiplying B by a diagonal 
matrix with 1 on the diagonal in the first 2N rows and 2 in the last N rows we have a W-matrix. Using 
the W-matrix Lemma 3.1 it therefore follows that small perturbations around an equilibrium solution are 
stable. 

Equivalent reaction rates can be defined as follows, 


k — 

^eq — 


N 

Wy fe U 

c 

i^j^k—1 


Ip.n — 


N 

^ ^ Lijk 


(3.52) 


Insert u(a;, t), v(a;, t), and w(a;,t) = w{x,t) Woc into (3.43)-(3.45) as in (3.38) and let 7 ^, = e’^JDwoo. As 
expected we recover the reaction-diffusion PDE for u, v, and w 


du 

'di 

dv 

dt 

dw 

'm 


7 u Au - keqUV + leg W, 
7 „ Av - keqUV + leg W, 
7 u, Aw + keqUV - leq W. 


(3.53) 

(3.54) 

(3.55) 


4. The internal states approximation of the FPDE. In this section we start off with the promising 
observation made in [38] that non-Markovian waiting times can be arbitrarily well approximated by a set of 
Markovian waiting times, each associated with its own internal state. In turn, those states are to be visited 
according to a certain random walk model which again can be taken as Markovian, all in all resulting in a 
computationally quite attractive modeling framework for subdiffusion and reactions. 

In Section 4.1, we recapitulate the basic internal states subdiffusive model and its relation to the FPDE, 
and in Section 4.2 we determine its asymptotic behavior for short and long times. In Section 4.3, we consider 
the same framework in the presence of reactions with two different mesoscopic models. In particular, we 
discuss the feasibility of obtaining coarse-grained macroscopic coefficients in a FPDE from observations of 
subdiffusive systems. 

4.1. Internal states diffusion system. The asymptotic behavior of the waiting time (2.9) at large 
time follows the power law 


V'(t) Ri Aa 



(4.1) 


with 


A 


a 


sin(7r a) 
TT 


r(l + a). 


(4.2) 


With a change of variable in the Euler’s P function, the diffusive representation of the totally monotone 
function in (4.1) is [7, 17, 20, 48] 


1 

ti+“ 


1 

r(l-ba) 


pOO 

/ e-" 

Jo 


ds. 


(4.3) 


The diffusive representation of (4.3) is approximated by using a quadrature formula in N points, with 
weights jli and abscissae sf. 


1 

ti+“ 


N 




t 


2=1 
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(4.4) 







leading to the diffusive approximation [2]. 

Our objective is then to approximate the function Fex{t) = by Fapprox{t) 
interval [tminffmax]- One possibility to quantify the error of the model £mod is 


N 

^ Pi in the time 

i=l 


^mod 


F„, 




Fexit) 






Fexit) 



(4.5) 


Based on the error (4.5), a nonlinear optimization is shown in [2] to be a better way to determine the 
parameters and Si than Gaussian quadrature. Consequently, this method is used in all what follows. 
Setting 


TV - \ 
Pi 


N 


i=l 


Si 


Si 




Si 


the approximation of the waiting time (4.1) is 


N 


^{t) Ri X! ^ ® 


(4.6) 


(4.7) 


The jump length variance T? (2.3) is then computed using (2.12). 

The waiting time PDF (2.9) and its asymptotic expansion (4.1) are compared in Figure 4.1. The 
parameters are those used in the numerical experiments in Section 5.2 (Table 5.1, second set of parameters). 
Anomalous diffusion is expected in the time range of interest [tmm, ^max]- In this time interval, the expansion 
(4.1) is already accurate as illustrated in Figure 4.1. 



Fig. 4.1. Section 4-1- Waiting time PDF (2.9) (blue solid line), and its asymptotic expansion (4.1) (red dotted line). The 
scales are logarithmic on both axes. 

The CTRW algorithm in Section 2 is extended to a multistate CTRW (MCTRW) algorithm in [38] for a 
diffusive system with internal states. One diffusive jump in the algorithm is performed as follows. Firstly, a 
state i is drawn with probability pi. Secondly, the waiting time is sampled from an exponential distribution 
with the PDF ijjiit) = The PDF of the jump length in an internal state i is \i{x) as in (2.4) with 

variance cr^. Finally, the length of the jump is drawn with the normal distribution of Xi. 
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The combined PDFs for the jump length and the waiting time of the internal states are 


N 

A(a;) = Kix), 


N N 

2=1 2=1 


(4.8) 


The reaction-diffusion system (3.13) in [38] is derived from the MCTRW algorithm and the approximation 
(4.7). In the *th internal state, the diffusion of Ui is ordinary, whereas the sum of the concentrations of the 
internal states at the macroscopic level t/ = u diffuses anomalously within the time range [tmimtmax]- 
A collection of internal states with waiting times ipi (t) are here approximated by one state with the waiting 
time 'ip{t) in (4.7). 


4.2. Diffusive behavior of the internal states model. To mimic a subdiffusive behavior and to 
define a Markov process, the waiting time PDF is approximated by a sum of N exponentials (4.7). Let 
us assume that ti < T 2 < ■ ■ ■ < T]\j. Equation (4.7) implies the following. 

• At small times, the Taylor expansion of '0(t) is 


V'(t) ~ vi(t) = 


N 

E 

2 = 1 


/i* A ^ e 


N 

E 


.-1 




1 - 




t^o+ 


(4.9) 


with Teq = ( X) ) / ( X) A” ) ■ At Small time, tp{t) is therefore equivalent to a Poisson law 

with parameter Teq. Hence, ordinary diffusion is expected for t <C ti; 

• As we see in (3.18), the macroscopic U satisfies a diffusion equation with a diffusion coefficient 7 (x, t) 
(3.19) varying in space and time. For large t, the long-time diffusion coefficient is given by (3.20). 
It does not depend on space and time. Consequently, a return to ordinary diffusion is expected. 

• For t G [Ti,TAr], a subdiffusive behavior is expected with subdiffusive exponent a. 

This is in agreement with observations in physical experiments [23, 25, 30, 41]. Numerical illustrations are 
found in Section 5. 


4.3. Internal states reaction-diffusion system. In this section, the mesoscopic model with N in¬ 
ternal states of the participating molecules in Section 3 is compared to the macroscopic FPDE models I and 
II in Section 2.3. The annihilation reaction, monomolecular reactions and bimolecular reactions summarized 
in Table 3.1 are investigated. For mesoscopic models with certain reactions, there are corresponding FPDE 
models but in other cases the macroscopic level with summation over the internal states is not so easily 
expressed as a FPDE. 

4.3.1. Annihilation process. We consider one species A with N different internal states. Each inter¬ 
nal state can be annihilated with rate ki 


A,-^0, i = (4.10) 

The internal states diffusion system (3.22) without v is then modified as follows 

du 

— = D Au-f A u — iTi u. (4-11) 

Let ki = k/Ti in (4.10) as in Table 3.1. Then Ki = kT and (4.11) can be written 

/9ll 

= T (cr^ Au — fc u) -I- A u. (4-12) 

at 

In model I (2.17), the fractional order temporal derivative g^i-c acts on both the Laplace operator and the 
reaction term. In the internal states diffusion system (4.11) and (4.12), the diffusion matrix is scaled by the 
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waiting time matrix T to mimic an anomalous behavior. Scaling the reaction matrix Ki in the same way 
will approximate the macroscopic model I. The waiting time for diffusive and reactive events is the 
same for all internal states. The relation between k in (4.12) and fc* in (2.17) is the same as between and 
Ka in (2.12), 


fc* (T^ 
~k ~ 1^ 


T . 


(4.13) 


For an approximation of model II, take K and Ki to be fc J in (3.26). This corresponds to an annihilation 
case in Table 3.1 with all = I. The right hand side of (4.11) is with this Ki 

Fill 

— = cr^ T Au + A u — fc u. (4-14) 

at 

The waiting time ipii't) for diffusive and reactive events is different for the internal states. 

If Ki and T and Ki and A commute then the equation for u in (4.14) can be written 


— = cr^ T Au + A u - fc u. (4.15) 

This equation corresponds to model II in Section 2.3.1 at the macroscopic level. The relation between fc in 
(4.14) and fc* in (2.24) is 


fc* = fc. (4-16) 

After a change of variables U = U in the model II FPDE (2.24), the FPDE (2.30) is obtained. In the 

same manner, introduce a change of variables u = u in (4.11). Then 

^ =a'^TAu + Au. (4.17) 

A sufficient condition for Ki to commute with T and A is that Ki = kl as is the case in (4.14). If Ki 
does not commute with T and A, the macroscopic equation for U = u may not be as simple as (2.30). 
Consider u, the total amount of A in the different states in 11, defined by 


u{t) = / u(a:, t) dn. 

Jq 

Integrating (4.11) and using the Neumann bonndary condition in (2.19) and (2.20) leads to 

—— = A u — Kl u. 

dt 

The time evolution of u is then 

u(t) = e(^-^i^‘u(0). 


(4.18) 


(4.19) 


(4.20) 


By Gerschgorin’s theorem for the eigenvalues of a matrix, the eigenvalues of A — Ki are strictly in the left 
half plane. By (4.19), the equation for the sum over all internal states U = e^u is 


with 


dt 


-e^ Kl u = -fc'(t) C/, 


k'{t) 


K\ u 


Kl ■^'^)*u(0) 
qT g(A-lCi) i 


(4.21) 
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(4.22) 




In model I with Ki = kT, k! varies in time depending on u, and the kinetics is then anomalous. Let Ai be 
the eigenvalue oi A — Ki with maximum real part and Si the corresponding eigenvector. Then for large t 


k'{t) 




(4.23) 


cf. (3.20). On the contrary, in model 11 k' = k does not depend on time, and the kinetics is then ordinary. 
This conclusion agrees with the comments in Section 2.3.1. 

4.3.2. Monomolecular reactions. We study two species A and B with N different internal states 
and the reversible reactions in (3.21). The internal states diffusion system is found in (3.22) and (3.23). The 
aim of this section is to choose the reaction matrices K and L in (3.26) and (3.27) to mimic either model 1 
(2.27) or model 11 (2.30) and to discuss other alternatives. 

The reaction rates are first chosen to be scaled by the waiting time matrix and one state i is transformed 
to the same state i in 


\z j 'T' 

Ai ^ Bi, i = l,...,N. 

tjri 


(4.24) 


Then Ka = fc/r^ and La = ijTi, cf. Table 3.1. With K = kT and L = iT in (3.22) and (3.23) we have 


— = T (tr^ Au —fcu + ^v) + Au, 

C/ L 

r\ 

= T (cr^ Av + fcu — £v) + Av. 


(4.25) 


The diffusion and the reactions have the same waiting time as in model 1 in (2.27). 

The steady state of model 1 in (4.25) has an analytical solution. Let Uoo span the nullspace of A (see 
(3.10)) and insert the constant solutions in space mUqo and vUoo with ||uoo||i = 1 as u and v in (4.25). 
Then the right hand sides are 


-kuTuoo + ivTuoo = 0 , , . 

kuTu^-ivTu^ =0. 

Both equations in (4.26) are satisfied ii ku = £v. Hence, the steady state solution is m(uoo , {k/£) Uqo), where 
u depends on the initial data. By mass conservation in Lemma 3.1 

e’^u(O) + e^v(O) = e^u(t) + e’^v(t) = it(e^Uoo + (fc/£)e^Uoo) = u{l + k/i). (4.27) 

The macroscopic reaction coefficients can be regarded as time and space dependent and the macroscopic 
equations can be expressed without fractional derivatives. Introduce small perturbations (5u and (5v around 
the steady state in (4.25) 


u(a;, t) = u Uoo + (5u(a;, t), 'v{x, t) = v Uqo + S'v{x, t). 

Then for u = m + (!l(||(5u||) and V = v = u + (!l(||(5v||) we derive 

dU 

T A(5u — kue^ T\ioo — k T 5\i + £v T Uoo + £e^ T 5^ = —K' U + L' V + (!l(||5u||), 
at 

dV 

T A^v + kue^ Tu^o +ke^ T 5u. — Ive^ T\iao — i T dv = K' U — L' V + C>(||(5v|l), 
at 

(4.28) 
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where 


e'^Tjuu^ + Su) 
^{uu^+Su) 
e^T{vUoo + Sv) y 
^{vu^ + S^r) 


^Tuoo + 0(||HI), 

Tuoo + O(||,5v||). 


(4.29) 


The effective macroscopic reaction coefficients are K' and L' when t ^ 0, cf. (3.36) with Ki = kT. 

With more general K and L in (3.21), the structure of the equations is the same as in (4.25) but there 
is no corresponding macroscopic FPDE. For example, let K = T ef^ T and L — Teg^ T with T e > 0 
and T e > 0. Then Ki = kT with k = T e and Li = iT with i = ^ T e. The equations for u and 
V are 


— = T (cr^ Au — fc u + g T v) + A u, 
at 

/9v 

— = T(cr2 Av + fe^Tu-£v) + Av. 


(4.30) 


The macroscopic waiting time for the reactions is not obvious in this case. 

As in (4.24), let one state i be transformed to the same state i and choose the reaction rates to be 
constant for all states Ka = k and La = t with = 1 in Table 3.1. Insert K = kl and L = il into (3.22) 
and (3.23) to obtain 


5u 

dt 

'm 


T Au + Au — fcu + ^v, 
cr^ T Av + Av + A:u — £v. 


(4.31) 


Only the diffusion has different waiting times in different states while the waiting times for the reactions are 
the same in all states. Let Mr and Md be the reaction matrix and the operator for the diffusion and the 
change of state in (4.31) 


Mr 


( -kl 

il \ 


^ (j^TA + A 

0 \ 

1 kl 

-u) 

, Md = 1 

. 0 

ct^tA + A j 


(4.32) 


There is a transformation from (4.14) via (4.15) to (4.17). Since Mr and Md commute in (4.32), there is a 
similar transformation for (4.31). As in (4.14), e^u and e^v approximate model II in (2.30). 

If the diffusion rate is independent of the state but the waiting time for the reactions depends on the 
state in (4.24) then Da = <J^ in (3.13), Ka = k/ri and La = ^Iri. The mesoscopic model is 


— = T {—k u + £v) + Au + cr^ Au, 
f)v 

— = T {k u — iv) + Av + Av. 


(4.33) 


We have the ordinary diffusion but the reactions behave anomalously at the macroscopic level with the 
waiting time depending on the state. 

4.3.3. Bimolecular reactions. We consider three species A, B and C in N different internal states 
and the chemical reactions 


Ai + Bi^Ci, * — 1, • 


,iV. 


Ai 
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(4.34) 




This is a less general set of reactions than in (3.42) in that only molecules in the same internal state react 
with each other. The internal states reaction-diffusion system (3.43)-(3.45) is then modified as follows for 
state i 


dui 

dt 


dv^ 

dt 


dwj 

dt 


N 

JDi ^Ui ^ ^ -^ij ^iii '^i 4 ” ; 

N 

D-i “1“ ^ ^ ^iii 4“ ^iii ; 

i=i 

N 

Di Awi + ^ Aij Wj + Km Ui Vi - Lm Wi. 

i=i 


The reaction coefficients in (4.35) are chosen to be Km = k/ri and Lm = ^iTi. Then at state i 


duj 

dt 


dvj 

dt 


dwi 

~dt 


1 


Ti 


1 


n 


1 


n 


N 

(cr^ Aui — kuiVi + Iwi) + ^ Aij Uj , 

i=i 

N 

(a^ Avi - kui Vi +£wi) +'^AijVj, 

j=i 

N 

(a^ Awi + kui Vi — £wi) + ^ Aij Wj. 

i=i 


(4.35) 


(4.36) 


The reactions and the diffusion have the same waiting time in each state in (4.36) as in (4.25). Hence, the 
equations will approximate model I at the macroscopic level in (2.31). By letting Km = k and Lm = i 
independent of the internal state in (4.35), equations similar to the model II equations (4.31) are obtained 
but there is no transformation (4.15) with a constant Ki. 

We have found that certain annihilation, monomolecular, and bimolecular reactions at the mesoscopic 
level have a macroscopic counterpart as a FPDE. Furthermore, the long time behavior of the reactions tends 
to that of ordinary reactions without internal states. This will be confirmed in numerical examples in the 
next section. 


5. Numerical experiments. Since this paper focuses on the mesoscopic approximation of subdiffu¬ 
sion, we will investigate only the ID reaction-diffusion system. However, the strategy proposed here can be 
extended straightforwardly to 2D and 3D geometries. 

The integration of the internal states systems (3.13), (3.22)-(3.23) and (3.43)-(3.45) is detailed in Sec¬ 
tion 5.1. The general configuration of the numerical experiments is introduced in Section 5.2; physical 
parameters, numerical parameters of the discretization, and initial conditions. Then the numerical exper¬ 
iments are described. In Section 5.3, both the diffusive approximation introduced in Section 4.1 and the 
numerical method introduced in Section 5.1 are verified in comparisons. In Sections 5.4 and 5.5, the method 
is applied to both model I and model H for reactive systems. The differences between the two models are 
also illustrated. Finally examples of bimolecular reactions, for which no analytical solutions are available, 
are presented in Section 5.6. 

5.1. Numerical modeling. In order to integrate the internal states systems (3.13), (3.22)-(3.23) and 
(3.43)-(3.45), a uniform grid is introduced with mesh size h and time step At. The approximation of the 
exact solution u is denoted by u" at Xj and t". The Laplace operators involved in the internal states 
systems (3.13), (3.22)-(3.23) and (3.43)-(3.45) are discretized using second order centered finite differences 
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as in (3.12): 


A 1 

^ - 2 u, + u,_i). 


(5.1) 


The jump coefficients in (3.12) are Xji = Xj 2 = 1/ft.^ and Xj = 2//i^. The resulting system of ODEs in time 
can be written 


dt 


Fi (Uj_i, Uj, Uj + i) + Fni (Uj) , 


(5.2) 


where Fi contains the discrete Laplacian and the linear reaction terms and the Fnf, the nonlinear reaction 
terms. The system of ODEs (5.2) is discretized using the following finite difference scheme 


n+1 


At 


- u” 1 


,n+l 


n+1 

3 




(5.3) 


In the case of an annihilation process and a monomolecular reaction, there is no nonlinear reaction term. 
Hence, Fng is zero, and (5.3) reduces to the Crank-Nicholson scheme. It is second-order accurate in space 
and time and it is unconditionally stable. In the case of a bimolecular reaction, because of the nonlinearities, 
the method is first order accurate in time. 


5.2. Configurations. In order to demonstrate the ability of the present method to be applied to a wide 
range of problems, we numerically test in Section 5.3 two different sets of parameters, given in Table 5.1. 
Then only the first set of parameters is used in Sections 5.4, 5.5 and 5.6. In our code, a, Ka, tmin, tmax are 
input parameters, and cr^, t, iV, Tj, are ouput parameters. The quadrature coefficients and /ii in (4.7) 
are determined by nonlinear optimization [2] and the corresponding model error Smod (4-5) is also given in 
the table. 

The computational domain is D = [—1,1] in Figure 5.1 and [—10,10] in Figure 5.2, is discretized with 
Nx = 128 grid points. Neumann boundary conditions are used. As the initial condition, we use a Gaussian 
g{x), centered at point (0,0) and of variance + = 10“^, rather than a Dirac distribution to avoid spurious 
numerical artifacts. Moreover, each internal state is initialized by the weight g.i. 



Parameters 

Set 1 


Set 2 


Physical parameters 

a. 

0.5 


0.5 



Ka (m2s-“) 

0.04 


0.04 



i'min (s) 

10-4 


10-3 



tmax (s) 

5 ■ 10 

-2 

1 



(m2) 

3.49 ■ 

10-4 

7.18 • 

10-4 


T (s) 

7.62 ■ 

10-® 

3.22 • 

10-4 

Optimization 

N 

4 


5 



Tl (s) 

9.51 ■ 

10-® 

7.58 ■ 

10-4 


T 2 (s) 

5.40 ■ 

10-4 

3.55 • 

10-3 


T3 (s) 

3.09 ■ 

10-3 

1.66 • 

10-2 


T 4 (s) 

2.13 ■ 

10-2 

7.89 • 

10-2 


T5 (s) 

- 


4.78 • 

10-1 



4.96 ■ 

10-4 

3.23 • 

10-4 


/^2 

2.07- 

10-4 

1.48 • 

10-4 


M3 

8.80 ■ 

10-2 

6.84 • 

10-2 


M4 

4.42 ■ 

10-2 

3.24 • 

10-2 



- 


1.85 ■ 

10-2 


^mod 

5.25 ■ 

10-2 

2.92 ■ 

10-2 


Table 5.1 


Parameters used in the numerical experiments. 
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5.3. SubdifFusion. The aim of the first test is to check the accuracy of the numerical method above 
when no reaction occurs. The following initial conditions are used 

u{x,0) = ^lg{x). (5.4) 

Figure 5.1 compares the numerical solution ?7 = u obtained with the PDFs in the internal state diffusion 
system (3.13) with the analytical solution of the FPDE (2.16). Figure 5.1-(a) corresponds to the first set of 
parameters, given in Table 5.1, at time = 5 • 10“^ s and Figure 5.1-(b) corresponds to the second set of 
parameters at time t 2 = 1.5 • 10“^ s. Note in Table 5.1 that the macroscopic parameters a and Ka, involved 
in (2.11), are the same for both sets of parameters. Consequently, they approximate the same macroscopic 
FPDE, but in two different time intervals [tmin,tmax]- In both cases, the final time of the simulation in 
Figure 5.1 is inside this interval. Excellent agreement is found between the two solutions. 

Two errors should be mentioned here: the modeling error, defined as the difference between the internal 
states diffusion model and the FPDE model, and the numerical error Snurm resulting from the numerical 
discretization of the internal states diffusion system. The modeling error is dependent on Smod in (4.5). The 
error Smod is given in Table 5.1 and the total error Stot is measured in Figure 5.1 as the difference between the 
numerical solution obtained with the PDF model (red circles) with the analytical solution obtained with the 
FPDE model (black solid line). In Figure 5.1-(a) Stot — 2.33 • 10“^ and in Figure 5.1-(b) Stot — 3.65 • 10“^. 


(a) (b) 




Fig. 5.1. Section 5.3. Comparison between the numerical values (circle) and the analytical values (solid line) of the 
concentration U = u of A. (a): Set 1 of parameters at ti = 5 ■ 10“® s, (b): Set 2 of parameters at t 2 = 1.5 • 10“^ s. 

Figure 5.2 shows the mean square displacement divided by the time, which is constant for ordinary 
diffusion, see (2.8). According to the analysis in Section 4.2, ordinary diffusion is observed for I ^ ^ min Te 

and for t ^ max t/. For intermediate times, the diffusion is anomalous by construction, cf. Section 4.1. 

e=i,...,N 

The subdiffusive exponent a is measured by linear regression: a ~ 0.4997 in Figure 5.2-(a) and a ~ 0.5454 
in Figure 5.2-(b). In both cases the measured subdiffusive exponent is close to the theoretical one. 


5.4. Annihilation process. The purpose of the second test is to investigate the accuracy of the 
numerical method in the case of an annihilation process. Two different models are possible in Section 2.3.1 


Model I A, ^ 0, 

Model II A, 0. 


(5.5) 
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Fig. 5.2. Section 5.3. Mean square displacement as a function of time, (a): Parameters set 1, (b): Parameters set 2. 
The scales are logarithmic on both axes. 



First 

run 

Second 

run 


Model I 

Model II 

Model I 

Model II 

k 

0.1 

8.727- s-i 

5.729 ■ 10-3 

50 s-i 

A:* 

8.727 ■ s-“ 

8.727 ■ 10“"' s-i 

50 s-“ 

50 s-i 


Table 5.2 

Reaction rates of annihilation process. 


The reaction coefficient k used in the simulations and the corresponding macroscopic reaction rate fc* (4.13)- 
(4.16) are given in Table 5.2. The rate k has been chosen such that the two models have the same macroscopic 
reaction rates in (4.13) and (4.16). Note the order-of-magnitude difference between the value of k in model 
I and II, which is explained by the fact that the units of k and fc* are are different. The initial conditions 
are the same as in the previous test. Figure 5.3 compares the numerical solution U = u obtained with 
the internal state reaction-diffusion system (4.11) with the analytical solutions of the FPDEs in (2.22) and 
(2.26) at time ti = 10“^ s. For both models, there is excellent agreement between the two solutions. The 
analytical solution is also given in the case where no reaction occurs (fc = 0). For a small macroscopic 
reaction rate (first run, top of Figure 5.3), no reaction takes place yet in the case of model II. For large 
macroscopic reaction rate (second run, bottom of Figure 5.3), all particles have disappeared in model I. The 
differences between the two models are clearly illustrated in Figure 5.3. 

The time development of k' (4.22) is plotted in Figure 5.4. In model I (first run in Table 5.2), k' is almost 
constant at small t and at large t, meaning that the kinetics is ordinary there but at intermediate t k' varies 
in time implying that the kinetics is anomalous. This is illustrated in Figure 5.5, where the time evolution 
of the total amount of kl, i.e. U in (4.21), is found. Exponential decay characterizing ordinary kinetics 
is observed for t < tmin and t > tmax- The exponential parameters are measured by linear regression in 
Figure 5.5 resulting in fc'(O) « 28.97 and k'^ « 674.08 when t —>• oo. In both cases, the measured exponential 
parameter is close to the theoretical one in (4.22) (blue values in Figure 5.4). For tmin <t< tmax, U does not 
decrease exponentially. In model II (second run in Table 5.2), Figure 5.4 shows that k' = k does not depend 
on time and exponential decay is observed for all times in Figure 5.5. The measured exponential parameter 
k' ~ 50.00 is close to the theoretical one. These observations agree with the analysis in Sections 2.3.1 and 
4.3.1. 
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Model I 


Model II 


1 

0.8 

0.6 

0.4 

0.2 


O' 

-1 


1 

0.8 

0.6 

0.4 

0.2 


O' 

-1 


-0.5 


0.5 


-1 


-0.5 



analytical FPDE 

1 


analytical FPDE 


without reaction 



without reaction 


analytical FPDE 



analytical FPDE 

! ‘ 

with reaction 



with reaction 


numerical PDE’s 



numerical PDE’s 

;' 

with reaction 



with reaction 

; ' 

; ; 

0.6 



Z) 



0.4 



0.2 




O' 



0.5 



analytical FPDE 

1 


analytical FPDE 


without reaction 



without reaction 


analytical FPDE 


I; 

analytical FPDE 

! ‘ 

with reaction 



with reaction 


^ numerical PDE’s 



numerical PDE’s 

;' 

with reaction 


;' 

with reaction 

; ' 

; ; 

0.6 

I 1 


Z) 

'A' 

i \ 

0.4 

■n' 

1 \ 

0.2 

1 \ 

\ 






O' 



-0.5 


0.5 


-1 


-0.5 


0.5 


Fig. 5.3. Section 5,4 Comparison between the numerical values (circles) and the analytical values (black line) of the 
conecntration U = u of A at time t = 10“^ s. Top: First run, bottom: Second run. 



Model I 

Model II 

k 


1.719 

10^ 

15 s 

-1 

1 


3.437 

103 

30 s 

-1 

fc* 


15 s 

— q: 

15 s 

-1 



30 s 

— q: 

30 s 

-1 

'tJ'OO 

(molm 

8.698 ■ 

10-3 

8.698 ■ 

10-3 


(mol ) 

4.349 ■ 

10-3 

4.349 ■ 

10-3 

keq 

(s-i) 

1.047 

106 

15 


ieq 

(s-1) 

2.094 

10® 

30 



Table 5.3 

Monomolecular reversible reaction. Reaction rates and theoretical steady state. 


5.5. Monomolecular reaction. The numerical method is here applied to monomolecular reversible 
reactions. Consider the reactions 


Model I 

fc/Tj 

Ai ^ , 

IjTi 

Model II 

A ^ 

Ai^Bi^ 


(5.6) 


where the coefficients k and i used in the simulations are giveir in Table 5.3. As discussed in Section 4.3.2, 
these cases correspond to the macroscopic FPDEs (2.27) and (2.30) with macroscopic reaction rates fc* and 
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Model I 


Model II 



t 

n 

nin 

nax 





Fig. 5.4. Section 5.4 k' in (4.22) in terms of time. Left: Model I, right: Model II. The scales are logarithmic on both axes. 



Fig. 5.5. Section 5.4- Total amount of A (4.21) in terms of time. Upper: Model I, lower: Model II. The scale is 
logarithmic on the y axis. 


given in Table 5.3. The rates k and I have been chosen such that the two models have the same macroscopic 
reaction rates A:* and and the same steady states. Figure 5.6 then illustrates the differences between the 
two models. The following initial conditions are used 

u(a;,0) = v(a;,0) =/X 5 (x). (5.7) 

Figure 5.6 compares the numerical solutions of C/ = u and F = e^v obtained with the internal state 
reaction-diffusion system (3.22)-(3.23) with the analytical solutions of the FPDEs in Section 2.3.2. The 
difference between the models is obvious in Figure 5.6-(a) corresponding to model I at time ti = 10“^ s and 
Figure 5.6-(b) with the model II at the same time. For both models, good agreement is found between the 
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mesoscopic PDE and the macroscopic FPDE solutions. 

The numerical solution of the internal states reaction-diffusion system at time t 2 = 10 s is depicted in 
Figure 5.7 when the steady state is reached. We note that the numerical values of the steady states are close 
to the theoretical ones, given by the kernel of the matrix B in (3.28). Using (3.36), the equivalent reactions 
rates keq and ^eq are computed, cf Table 5.3, and the property k^q Uoo = ^eqVoo when t —>■ oo in (3.40) and 
(3.41) is verified. 


Model I 


Model II 



-1 -0.5 0 0.5 1 -1 -0.5 0 0.5 1 

X X 


Fig. 5.6. Section 5.5. Comparison between the numerical values (cireles and triangles) and the analytical values (solid 
line and dashed line) of the concentration U = u, F = v of A and B at time t\ = 10“^ s. 


Model I 


-0.5 


Model II 


0.5 


-0.5 


X 10'^ 



X 10'^ 






lO"^ 


8.698- 

lO'^ 



- U(x,t) 

- V(x,t) 
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7 

> 

zS 
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- U(x,t) 
V(x,t) 





_ 


5 

4.349- 

lO’"® 






0.5 


Fig. 5.7. Section 5.5. Numerical values of the concentration U = u, F = v of A and B at time t 2 = 10 s. 


5.6. Bimolecular reaction. The purpose of the last example is to establish whether the numerical 
methods presented in this paper can be used to handle more complex reactions. As an example, we consider 
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Reaction I 

Reaction II 

Reaction III 

k 

37.5 mmol“^ 

37.5 mmol-^ 

7500 mmol-3 

e 

0.25 

0.25 

5 s-3 

Uoo (molm“^) 

8.204- 10-3 

5.826 ■ 10-3 

3.772 ■ 10-3 

Voo (molm“l) 

6.573 ■ 10-3 

4.195 ■ 10-3 

2.141 ■ 10-3 

Woo (molm“l) 

1.582 ■ 10-3 

3.959 ■ 10-3 

6.014- 10-3 

keq (mmol~^s~^) 

2.450 ■ 103 

1.187 ■ 10® 

2.483 - IQi 

£eq (s 

8.351 ■ 103 

7.328 ■ 10^ 

10^ 


Table 5.4 

Bimolecular reversible reaction. Reaction rates and theoretical steady states. 


the following bimolecular reversible reactions 



Reaction I 

kjri 

Ai ~\- Bi ^ Ci, 
eir, 


Reaction II 

Ai + Bj ,— Ck, 

t/Tk 

(5.8) 

Reaction III 

Ai Bi^Ci. 



I 


The reaction rates are found in Table 3.1 with for reaction II chosen as in (3.7) with 0 = 1/2 since A and 
B diffuse with the same speed. The rates of reactions I and III are as in model I and II in Sections 4.3.1 
and 4.3.2 but reaction II is more general. The initial conditions are 

u(a;,0) = i/xg(a;), v{x,0) = w(x, 0) =/xg(a;). (5.9) 

The numerical solutions U = u, V = v and IT = w of the internal states reaction-diffusion system 
are shown in Figure 5.6 at time ti = 10“^ s (top) and at time ^2 = 10 s (bottom) when the steady state 
is reached. For reaction I, we expect that the internal states model in (4.36) approximates the macroscopic 
FPDE model I (2.31). No analytical solution of the FPDE is available in this case. For reaction II and 
reaction III, the macroscopic level with summation over the internal states is not so easily expressed as a 
FPDE. Using (3.48)-(3.52), the equivalent reactions rates keq and igq are computed in Table 5.4 and the 
property keq Uoo 'Coo = ■^eq Woo when f —oo in (3.53)-(3.55) is verified. 

6. Conclusion. A numerical method is presented here for simulating fractional-in-time reaction-diffusion 
equations. A diffusive representation transforms the function involved in the mesoscopic CTRW on 

a lattice, into a continuum of decreasing exponentials, approximated by quadrature formulae. The CTRW 
model is then replaced by an approximation, much more tractable numerically. Contrary to the approach 
used in [38], the coefficients of the diffusive approximation are determined by a nonlinear optimization pro¬ 
cedure, leading to a smaller number of internal states. At the macroscopic level, the internal states diffusion 
system thus obtained corresponds to the fractional diffusion equation in a chosen time interval. In contrast 
to the EPDE model, the diffusion in the internal states model is ordinary at small and large times, but it 
is anomalous at intermediate times. This behavior can also be observed in crowded systems of hard-spheres 
due to caging effects, and hence the model used herein may actually be a better model for microscopic 
crowding than the traditional EPDE description. 

The internal states model for diffusion in [38] is here extended to account for chemical reactions. On the 
macroscopic FPDE level, two different models for reactions with subdiffusion are investigated. In model I 
the fractional derivative acts on both on the standard diffusion term and the reaction term, whereas in model 
II the fractional derivative acts only on the diffusion term. Both macroscopic FPDE models correspond to 
a mesoscopic internal states model with particular reaction coefficients. However, the opposite is not true; 
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Reaction I 


Reaction II 


reaction III 



XXX 


Fig. 5.8. Section 5.6. Numerical values of the concentration U = u, V = e^v, W = w of A, B and C at time 
ti = 1 ■ 10“^ s (top) and = 10 s (bottom). 


mesoscopic models with general reactions may not have a simple interpretation at the macroscopic level. In 
model I, the reactions are subdiffusion controlled, that is the reaction kinetics is ordinary at small and large 
times, whereas it is anomalous at intermediate times. In model II, the reaction kinetics is ordinary for all 
times. Which one of these models provides a better description of a reaction system subject to subdiffusion 
does not have a simple answer. In either case, the present work provids a theoretical foundation for practical 
and efficient mesoscopic simulations. 
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Appendix A. Special functions. 

The Eox-H function is defined as a Mellin-Barnes integral [12] 


Tjm.n 



(oi, Ai) 

(02,^2) 

(ap,x 4 p) 



(blyBl) 

(62,52) 

{bq, Bg) 

Jl 


n T{bj - Bj s) X n r(l - aj + Aj s) 

j=i 

n r(l - bj + Bj s) X n r(aj - Aj s) 

j—rn-\-l j=n-\-l 


’’ ds. 


(A.l) 

where L is a certain contour separating the poles of the two factors in the numerator. The special case for 
which the Fox-H function reduces to the Meijer-G function is Aj = Bp = C, C > 0 for j = 1, • • • ,p and 
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Tjm.n 



(ai,C) 

(a2,C') •• 

• {ap, C) 

^ i 


ibi,C) 

ib2,C) •• 

■ NC) \ 

2j7rC A 


n T{bj - s) X n r(l - aj + s) 

Zzi_izi_ Z-/C 

n r(l-6j+s)x n r(aj-s) 

j=n+l 


~ C 


Mc 


Oj\ , • • • , dp 


(A.2) 
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